$NOFLOATCALLS
$NODEBUG
$STORAGE:2
c***********************************************************************

	subroutine rkp(x,ys,h,neqn,imeth,kt,est,yold,ynew,w,func)


	implicit double precision (a-h,k,p-z)
	external func
	dimension w(13,neqn),kt(neqn),ys(neqn),est(neqn)
	dimension al(12),b(12,12),yold(neqn),ynew(neqn)
	common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
     &		       a12,a13,er1,er3,er4,er5,er6,er7,er8,
     &		       er9,er10,er11,er12,er13,inum
	xs=x
	do 5 i=1,neqn
	   ys(i)=yold(i)
5	continue
	call func(xs,ys,kt,neqn)
	do 7 i=1,neqn
	   w(1,i)=kt(i)
7	continue

	do 20 j=2,inum
	   j1=j-1
	   xs=x+h*al(j1)
	   do 15 i=1,neqn
	      ksum=0.d0
	      do 10 m=1,j1
		 ksum=ksum+b(m,j1)*w(m,i)
10	      continue
	      ys(i)=yold(i)+h*ksum
15	   continue
	   call func(xs,ys,kt,neqn)
	   do 17 i=1,neqn
	      w(j,i)=kt(i)
17	   continue
20	continue

c****	evaluate YNEW[i] and the error estimates EST[i]

	if(imeth.EQ.1) then
	do 30 i=1,neqn
	   ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
     &					     +a5*w(5,i)+a6*w(6,i))
	   est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
     &				  +er6*w(6,i))
30	continue
	end if

	if(imeth.EQ.2) then
	do 40 i=1,neqn
	   ynew(i)=yold(i)+h*(a1*w(1,i)+a3*w(3,i)+a4*w(4,i)
     &					     +a5*w(5,i)+a7*w(7,i)
     &					     +a8*w(8,i))
	   est(i)=h*(er1*w(1,i)+er3*w(3,i)+er4*w(4,i)+er5*w(5,i)
     &				  +er6*w(6,i)+er7*w(7,i)+er8*w(8,i))
40	continue
	end if

	if(imeth.EQ.3) then
	do 50 i=1,neqn
	   ynew(i)=yold(i)+h*(a1*w(1,i)+a6*w(6,i)+a7*w(7,i)
     &			     +a8*w(8,i)+a9*w(9,i)
     &			     +a12*w(12,i)+a13*w(13,i))
	   est(i)=h*(er1*w(1,i)
     &		    +er6*w(6,i)+er7*w(7,i)+er8*w(8,i)+er9*w(9,i)
     &		    +er10*w(10,i)+er11*w(11,i)+er12*w(12,i)
     &		    +er13*w(13,i))
50	continue
	end if

	return
	end
